!=================Programa que calcula algumas derivadas===============

PROGRAM exerA

!======================================================================

IMPLICIT NONE

!========================Declaração de variáveis=======================

      INTEGER :: num_h, I, J, min_err
      REAL(8) :: f, der_f, sndder_f, trdder_f
      REAL(8), ALLOCATABLE, DIMENSION(:) :: h_values
      REAL(8), PARAMETER :: x0 = (1.0d0/3.0d0)
      REAL(8), ALLOCATABLE, DIMENSION(:,:) :: func_value
      !Define o formato para escrever na saída
      CHARACTER(LEN=*), PARAMETER :: FMT1 = "(7ES14.7)"

!=====================Abertura e leitura de arquivos===================

      OPEN(UNIT=10, FILE='tab1_in.dat', STATUS='old')
      OPEN(UNIT=20, FILE='tab1_out.dat', STATUS='unknown')
      
      !Lê o número de valores de h
      READ(10,*) num_h

      !Aloca a memória do vetor que terá os valores de h e da matriz
      !que é a tabela a ser escrita no arquivo de saída
      ALLOCATE(h_values(num_h),func_value(7,num_h))

      !Lê do arquivo para o vetor
      READ(10,*) h_values

      func_value = 0

!=========Calcula as derivadas para todos os valores de h==============

      DO I=1,num_h
            !Escreve o h atual na linha da matriz
            func_value(1,i) = h_values(I)
            !Simétrica de 3 pontos
            func_value(2,i) = ABS((f(x0+h_values(I)) - f(x0-h_values(I)))&
            &/ (2*h_values(I)) - der_f(x0))
            !Para frente de 2 pontos
            func_value(3,i) = ABS((f(x0+h_values(I)) - f(x0))/h_values(I)&
            &- der_f(x0))
            !Para trás de 2 pontos
            func_value(4,i) = ABS((f(x0) - f(x0-h_values(I)))/h_values(I)&
            &- der_f(x0))
            !Simétrica de 5 pontos
            func_value(5,i) = ABS((f(x0-2*h_values(I))-8*f(x0-h_values(I))&
            &+8*f(x0+h_values(I))-f(x0+2*h_values(I)))/(12*h_values(I))&
            &-der_f(x0))
            !Segunda simétrica de 5 pontos
            func_value(6,i) = ABS((-f(x0+2*h_values(I))+16*f(x0+&
            &h_values(I))-30*f(x0)+16*f(x0-h_values(I))-f(x0-2*&
            &h_values(I)))/(12*h_values(I)**2) - sndder_f(x0))
            !Terceira anti-simétrica 5 pontos
            func_value(7,i) = ABS((f(x0+2*h_values(I))-2*f(x0+h_values(I))&
            &+2*f(x0-h_values(I))-f(x0-2*h_values(I)))/(2*h_values(I)**3)&
            &-trdder_f(x0))
      END DO

      !Escreve no arquivo a tabela
      WRITE(20,FMT1) func_value

!===============Encontra o melhor h percorrendo a matriz===============

      DO I=2,7
            DO J=1,(num_h-1)
                  IF (func_value(I,J+1) < func_value(I,J)) THEN
                        min_err = J+1
                  END IF
            END DO
            WRITE(*,*) "No", I-1, "º caso, o melhor h é:", &
            &func_value(1,min_err)
      END DO

!==========================Fecha os arquivos===========================

      CLOSE(UNIT=10)
      CLOSE(UNIT=20)

!======================================================================

END PROGRAM exerA

!======================================================================


!========================Subroutines e funções=========================

!Valor de f no ponto x
REAL(8) FUNCTION f(x)
      IMPLICIT NONE

      REAL(8) :: x

      f = DEXP(4.0d0*x)*DCOS(x/2.0d0)
END FUNCTION

!Valor da primeira derivada de f no ponto x
REAL(8) FUNCTION der_f(x)
      IMPLICIT NONE

      REAL(8) :: x

      der_f = 4*DEXP(4.0d0*x)*DCOS(x/2.0d0)-(1.0d0/2.0d0)*DEXP(4.0d0*x)&
      &*DSIN(x/2.0d0)
END FUNCTION der_f

!Valor da segunda derivada de f no ponto x
REAL(8) FUNCTION sndder_f(x)
      IMPLICIT NONE

      REAL(8) :: x

      sndder_f = (63.0d0/4.0d0)*DEXP(4.0d0*x)*DCOS(x/2.0d0)-4.0d0&
      &*DEXP(4.0d0*x)*DSIN(x/2.0d0)
END FUNCTION sndder_f

!Valor da terceira derivada de f no ponto x
REAL(8) FUNCTION trdder_f(x)
      IMPLICIT NONE

      REAL(8) :: x

      trdder_f = 61.0d0*DEXP(4.0d0*x)*DCOS(x/2.0d0)-(191.0d0/8.0d0)&
      &*DEXP(4.0d0*x)*DSIN(x/2.0d0)
END FUNCTION trdder_f
